
%% This File Runs the Constant-Parameter VAR models %%

% When using this code or parts of it, please cite Paul, Pascal. 2019. 
% "The Time-Varying Effect of Monetary Policy on Asset Prices". 
% Review of Economics & Statistics.

clear all
clc

%% Choose Figure to Run

choose_figure       = 2;
colored_figures     = 1;
    % Figure 2:  VARX
    % Appendix A.14.
    % Figure 8:  VARX Wild Bootstrap
    % Figure 9:  Cholesky Decomposition
    % Figure 10: VARX excluding ELB
    % Figure 11: VARX Excess Bond Premium
    % Figure 12: VARX short-run restriction CPI
    % Figure 13: External instrument approach

figure_settings

% VAR Settings
p      = 4;         % Number of lags
T      = 400;       % Computed horizon IRFs
IRF_T  = 40;        % Plotted horizon IRFs
M      = 1000;      % Number of bootstrap repetitions
L      = 20;        % Block-size MBB

%% Load Data

cd '..\Data\Macro Time Series'

% Load Macro series (1960 M1 - 2017 M9)
macro           = xlsread('Macro Data.xls','Data','C2:I694');

% Convert stock prices, dividends, and house prices from nominal to real values
macro(:,2)      = 100*macro(:,2)./macro(:,5);
macro(:,3)      = 100*macro(:,3)./macro(:,5);
macro(:,4)      = 100*macro(:,4)./macro(:,5);

cd '..\Monetary Policy Surprises'

% Load Monetary Policy Surprises
surprises   = xlsread('Monetary_Surprises.xls','Scheduled_Meetings','C2:C405');

cd '..\..\Constant Parameter VAR'

%% Define VAR

vars            = [macro(:,1), macro(:,2), macro(:,3), macro(:,4), macro(:,5), macro(:,6)];
labels          = {'Federal Funds Rate','Stock Prices','Dividends','House Prices','Consumer Price Index','Industrial Production'};
taking_logs     = [0, 1, 1, 1, 1, 1];
differencing    = [0, 1, 1, 1, 1, 1];
if ebp == 1
    vars            = [macro(:,1), macro(:,2), macro(:,3), macro(:,4), macro(:,5), macro(:,6), macro(:,7)];
    labels          = {'Federal Funds Rate','Stock Prices','Dividends','House Prices','Consumer Price Index','Industrial Production','Excess Bond Premium'};
    taking_logs     = [0, 1, 1, 1, 1, 1, 0];
    differencing    = [0, 1, 1, 1, 1, 1, 0]; 
end
n_vars          = size(vars,2);

start_sample    = 12*(start_year-1960)+start_month;
end_sample      = size(vars,1)-12*(2017-end_year)-(9-end_month);
vars_y          = vars(start_sample:end_sample,:);

%% Pre-processing

% Pre-processing VAR
shorten_sample = 0;
for jj = 1:size(vars_y,2)
    if taking_logs(1,jj) == 1
        vars_y(:,jj) = 100*log(vars_y(:,jj));
    end
    if differencing(1,jj) == 1
        vars_y(:,jj)    = [0;diff(vars_y(:,jj))];
        shorten_sample  = 1 ;
    end
end

if shorten_sample == 1
    vars_y          = vars_y(2:end,:);
    start_sample    = start_sample + 1;
end

% Instrument on chosen sample
start_instr    = 12*(start_year-1984) + (start_month-1) + p + 1;
end_instr      = size(surprises,1)-12*(2017-end_year)-(9-end_month);
instr          = surprises(start_instr:end_instr,1);

%% Estimate VARX

Y = vars_y(p+1:end,:);
X = lagmatrix(vars_y,1:p);
X = X(p+1:end,:);

% Orthogonalize instrument to lags of Y
beta_instr  = X\instr;
instr       = instr-X*beta_instr;

if zero_restriction == 1
    bet_1 = [X instr ones(length(X),1)]\Y(:,1);
    bet_2 = [X instr ones(length(X),1)]\Y(:,2);
    bet_3 = [X instr ones(length(X),1)]\Y(:,3);
    bet_4 = [X instr ones(length(X),1)]\Y(:,4);
    bet_5 = [X ones(length(X),1)]\Y(:,5);
    bet_5 = [bet_5(1:end-1,1);0;bet_5(end,1)];
    bet_6 = [X instr ones(length(X),1)]\Y(:,6);
    bet   = [bet_1 bet_2 bet_3 bet_4 bet_5 bet_6];
else
    bet   = [X instr ones(length(X),1)]\Y;
end

% Generate residuals
res       = Y-[X instr ones(length(X),1)]*bet;

%% Estimate VAR with External Instrument Approach

bet_ident           = [X ones(length(X),1)]\Y;
res_ident           = Y-[X ones(length(X),1)]*bet_ident;
[S,F_ext_instr]     = ident_ext_instr(n_vars, p, instr, res_ident);
norm_shock          = S(1,1)/bet(end-1,1);

if zero_restriction == 1
    S = [S(1:4,1);0;S(6,1)];
end

%% Compute Impulse Responses

irfs      = irf_exog_var(bet, Y, T, p, res, norm_shock);
irfs_gk   = irf_ext_instr(bet_ident, Y, T, S, p);

%% Compute Confidence Bands via Wild Bootstrap or Moving-Block Bootstrap

% Generate disturbances for residuals and proxy
if Bootstrap == 1 % Recursive Wild Bootstrap
    eb = 1 - 2*(rand(size(res,1),M)>0.5);
elseif Bootstrap == 2 % Moving-Block Bootstrap
    [Res_m,Res_m_gk,Instr_m] = MBB(L,res,res_ident,instr,M);
end

% Store the IRFs
irf_store       = zeros(size(Y,2),T,M);
irf_store_gk    = zeros(size(Y,2),T,M);

for mm = 1:M
                
    if Bootstrap == 1
        eb_m    = eb(:,mm);
        res_m   = res.*(eb_m*ones(1,size(res,2)));
        instr_m = instr.*(eb_m*ones(1,size(instr,2)));
    elseif Bootstrap == 2
        res_m       = Res_m(:,:,mm);
        instr_m     = Instr_m(:,:,mm);
    end
    
    % Generate new series recursively with OLS estimates and bootstrapped residuals 
    % Set the first p observations to equal the data for simulation purposes
    gens        = zeros(size(res_m,1)+p,size(res_m,2));       
    gens(1:p,:) = vars_y(1:p,:);
    
    for tt = 1:size(res_m,1)
        xx = gens(tt:tt+p-1,:);
        xx = flipud(xx);
        xx = xx';
        xx = reshape(xx,1,size(xx,1)*size(xx,2));
        gens(p+tt,:) = [xx instr_m(tt,1) 1]*bet + res_m(tt,:);
    end
    
    % Use the generated data to re-estimate the VAR
    Y_m = gens(p+1:end,:);
    X_m = lagmatrix(gens,1:p);
    X_m = X_m(p+1:end,:);

    % Clean off the instrument
    beta_instr_m    = X_m\instr_m;
    instr_m         = instr_m-X_m*beta_instr_m;
    
    if zero_restriction == 1
        bet_1_m = [X_m instr_m ones(length(X_m),1)]\Y_m(:,1);
        bet_2_m = [X_m instr_m ones(length(X_m),1)]\Y_m(:,2);
        bet_3_m = [X_m instr_m ones(length(X_m),1)]\Y_m(:,3);
        bet_4_m = [X_m instr_m ones(length(X_m),1)]\Y_m(:,4);
        bet_5_m = [X_m ones(length(X_m),1)]\Y_m(:,5);
        bet_6_m = [X_m instr_m ones(length(X_m),1)]\Y_m(:,6);
        bet_5_m = [bet_5_m(1:end-1,1);0;bet_5_m(end,1)];
        bet_m   = [bet_1_m bet_2_m bet_3_m bet_4_m bet_5_m bet_6_m];
    else
        bet_m     = [X_m instr_m ones(length(X_m),1)]\Y_m;
    end
    
    res_new = Y_m - [X_m instr_m ones(length(X_m),1)]*bet_m;
    
    % Recompute VAR for external instrument approach
    if Bootstrap == 1
        res_m_gk       = res_ident.*(eb_m*ones(1,size(res_ident,2)));
    elseif Bootstrap == 2
        res_m_gk       = Res_m_gk(:,:,mm);
        instr_m_gk     = Instr_m(:,:,mm);
    end
    
    gens_gk         = zeros(size(res_m_gk,1)+p, size(res_m_gk,2));
    gens_gk(1:p,:)  = vars_y(1:p,:);
    
    for tt = 1:size(res_m,1)
        xx = gens_gk(tt:tt+p-1,:);
        xx = flipud(xx);
        xx = xx';
        xx = reshape(xx,1,size(xx,1)*size(xx,2));
        gens_gk(p+tt,:) = [xx 1]*bet_ident + res_m_gk(tt,:);
    end
    
    Y_m_gk = gens_gk(p+1:end,:);
    X_m_gk = lagmatrix(gens_gk,1:p);
    X_m_gk = X_m_gk(p+1:end,:);
    
    bet_ident_m   = [X_m_gk ones(length(X_m_gk),1)]\Y_m_gk;
    res_ident_m   = Y_m_gk-[X_m_gk ones(length(X_m_gk),1)]*bet_ident_m;
    [S_m,~]       = ident_ext_instr(n_vars, p, instr_m, res_ident_m);
    
    if zero_restriction == 1
        S_m = [S_m(1:4,1);0;S_m(6,1)];
    end
    
    norm_shock_m = S_m(1,1)/bet_m(end-1,1);
    
    % Recompute IRF
    irf_store(:,:,mm)       = irf_exog_var(bet_m, Y_m, T, p, res_new, norm_shock_m);
    irf_store_gk(:,:,mm)    = irf_ext_instr(bet_ident_m, Y_m_gk, T, S_m, p);
    
end

% Cumulating IRFs for first-differenced variables

for vv=1:size(vars_y,2)
    if differencing(1,vv)==1
        irf_store(vv,:,:)           = cumsum(irf_store(vv,:,:),2);
        irf_store_gk(vv,:,:)        = cumsum(irf_store_gk(vv,:,:),2);
        irfs(vv,:)                  = cumsum(irfs(vv,:),2);
        irfs_gk(vv,:)               = cumsum(irfs_gk(vv,:),2);
    end
end
            
% Take percentiles
            
irf_lb(:,:)  = prctile(irf_store(:,:,:),2.5,3);
irf_lbb(:,:) = prctile(irf_store(:,:,:),16,3);
irf_ubb(:,:) = prctile(irf_store(:,:,:),84,3);
irf_ub(:,:)  = prctile(irf_store(:,:,:),97.5,3);

irf_lb_gk(:,:)  = prctile(irf_store_gk(:,:,:),2.5,3);
irf_lbb_gk(:,:) = prctile(irf_store_gk(:,:,:),16,3);
irf_ubb_gk(:,:) = prctile(irf_store_gk(:,:,:),84,3);
irf_ub_gk(:,:)  = prctile(irf_store_gk(:,:,:),97.5,3);

irf_med(:,:)    = prctile(irf_store(:,:,:),50,3);
irf_med_gk(:,:) = prctile(irf_store_gk(:,:,:),50,3);

%% Compute Real Interest Rate and Fundamental Stock Price Response

% Computing IRF for real interest rate
irf_rr      = zeros(1,T);
irf_rr_lb   = zeros(1,T);
irf_rr_ub   = zeros(1,T);

irf_rr_gk      = zeros(1,T);
irf_rr_lb_gk   = zeros(1,T);
irf_rr_ub_gk   = zeros(1,T);

for jj=1:1:size(irf_rr,2)-12
    irf_rr(1,jj)    = irfs(1,jj)    - (irfs(5,jj+12)-irfs(5,jj)) ;
    irf_rr_lb(1,jj) = irf_lb(1,jj)  - (irf_lb(5,jj+12)-irf_lb(5,jj)) ;
    irf_rr_ub(1,jj) = irf_ub(1,jj)  - (irf_ub(5,jj+12)-irf_ub(5,jj)) ;

    irf_rr_gk(1,jj)    = irfs_gk(1,jj)    - (irfs_gk(5,jj+12)-irfs_gk(5,jj)) ;
    irf_rr_lb_gk(1,jj) = irf_lb_gk(1,jj)  - (irf_lb_gk(5,jj+12)-irf_lb_gk(5,jj)) ;
    irf_rr_ub_gk(1,jj) = irf_ub_gk(1,jj)  - (irf_ub_gk(5,jj+12)-irf_ub_gk(5,jj)) ;
end

irfs    = [irfs; irf_rr];
irf_lb  = [irf_lb; irf_rr_lb];
irf_ub  = [irf_ub; irf_rr_ub];

irfs_gk    = [irfs_gk; irf_rr_gk];
irf_lb_gk  = [irf_lb_gk; irf_rr_lb_gk];
irf_ub_gk  = [irf_ub_gk; irf_rr_ub_gk];

% Computing the fundamental stock price

indR    = n_vars + 1;
Lambda  = 1/(1+0.04/12);
qf      = zeros(size(irfs,2),1);
qf_gk   = zeros(size(irfs,2),1);
J       = T-150;

for k=0:IRF_T

     Funda      = zeros(J+1,1);
     Funda_gk   = zeros(J+1,1);

     for j=0:J
          Funda(j+1,1)      = Lambda^j*[(1-Lambda)*irfs(3,k+j+2)-irfs(7,k+j+1)/12];
          Funda_gk(j+1,1)   = Lambda^j*[(1-Lambda)*irfs_gk(3,k+j+2)-irfs_gk(7,k+j+1)/12];
     end

     qf(k+1,1)      = sum(Funda);
     qf_gk(k+1,1)   = sum(Funda_gk);

end

%% Bounds for Figures

figure_bounds

%% Figure: Exogenous Variable Approach

if choose_figure == 2 || choose_figure == 8 || choose_figure == 10 || choose_figure == 12

    figure(1)
    for mm=1:6
        subplot(3,2,mm)
        plot(irf_med(mm,:),'Linewidth',3,'Color',rgb(irf_color))
        hold on
        plot(irf_lb(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_ub(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_lbb(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_ubb(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        ciplot(irf_lb(mm,:),irf_ub(mm,:),band_color);
        ciplot(irf_lbb(mm,:),irf_ubb(mm,:),band_color);
        if mm == 1
            plot(irfs(n_vars+1,:),'Linewidth',3,'Linestyle','-.','Color',rgb(dash_color))
        elseif mm == 2
            plot(qf,'Linewidth',3,'Linestyle','-.','Color',rgb(dash_color))
        end
        plot(zeros(size(irfs(mm,:),2),1),'k')
        axis tight
        xlim([1 IRF_T])
        ylim([y_lim(mm,1) y_lim(mm,2)])
        t = title(labels(1,mm));
        set(gca,'Box','on','FontSize', 16)
        ylabel('Percent','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
        if mm == 5 || mm == 6
            xlabel('Months','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
        end
        set(t,'fontsize',24,'FontName','Palatino','Interpreter','Latex')
        grid on
        set(gca,'GridLineStyle','--')
        ax = gca;
        ax.GridAlpha = .5;     
    end
    h = gcf;
    set(h,'PaperPositionMode','auto')
    set(h,'PaperOrientation','landscape')
    set(h,'Position',[0 0 1500 1000])
    tightfig

    if choose_figure == 2
        print(h,'-depsc','..\Output\Figure_2')
    elseif choose_figure == 8
        print(h,'-depsc','..\Output\Figure_8')
    elseif choose_figure == 10
        print(h,'-depsc','..\Output\Figure_10')
    elseif choose_figure == 12
        print(h,'-depsc','..\Output\Figure_12')    
    end
end

%% Figure: External Instrument Approach

if choose_figure == 13

    figure(1)
    for mm=1:size(irf_med_gk,1)
        subplot(3,2,mm)
        plot(irf_med_gk(mm,:),'Linewidth',3,'Color',rgb(irf_color))
        hold on
        plot(irf_lb_gk(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_ub_gk(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_lbb_gk(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_ubb_gk(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        ciplot(irf_lb_gk(mm,:),irf_ub_gk(mm,:), band_color);
        ciplot(irf_lbb_gk(mm,:),irf_ubb_gk(mm,:), band_color);
        if mm == 1
            plot(irfs_gk(n_vars+1,:),'Linewidth',3,'Linestyle','-.','Color',rgb(dash_color))
        elseif mm == 2
            plot(qf_gk,'Linewidth',3,'Linestyle','-.','Color',rgb(dash_color))
        end
        plot(zeros(size(irfs(1,:),2),1),'k')
        axis tight
        xlim([1 IRF_T])
        ylim([y_lim(mm,1) y_lim(mm,2)])
        t=title(labels(1,mm));
        set(gca,'Box','on','FontSize', 16)
        ylabel('Percent','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
        if mm == 5 || mm == 6
            xlabel('Months','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
        end
        set(t,'fontsize',24,'FontName','Palatino','Interpreter','Latex')
        grid on
        set(gca,'GridLineStyle','--')
        ax = gca;
        ax.GridAlpha = .5;
    end
    h = gcf;
    set(h, 'PaperPositionMode','auto')
    set(h,'PaperOrientation','landscape')
    set(h,'Position',[0 0 1500 1000])
    tightfig
    print(h,'-depsc','..\Output\Figure_13')
end

%% Figure: Cholesky Decomposition

if choose_figure == 9
    
    cd 'Cholesky VAR'

    vars_chol           = [vars_y(:,6) vars_y(:,5) vars_y(:,3) vars_y(:,1) vars_y(:,2) vars_y(:,4)];
    vars_chol_labels    = {'Industrial Production','Consumer Price Index','Dividends','Federal Funds Rate','Stock Prices','House Prices'};
    differencing_chol   = [differencing(1,6), differencing(1,5), differencing(1,3), differencing(1,1), differencing(1,2), differencing(1,4)];
    index_Div           = 3;
    index_CPI           = 2;
    index_i             = 4;

    results             = VARmodel(vars_chol,p,1);

    [IRF,IRFopt]        = VARir(results,T,'oir',0);

    % irf_store_chol(t,k,m) where t=time, k=variable, m=sample
    [irf_store_chol]    = VARirband(results,IRFopt,M,95,'bs',index_i);

    % Cumulating IRFs for first-differenced variables over t
    for vv=1:size(vars_chol,2)
        if differencing_chol(1,vv) == 1
            irf_store_chol(:,vv,:)      =cumsum(irf_store_chol(:,vv,:),1);
        end
    end

    irfs_chol(:,:)       = prctile(irf_store_chol(:,:,:),50,3);
    irfs_lb_chol(:,:)    = prctile(irf_store_chol(:,:,:),2.5,3);
    irfs_ub_chol(:,:)    = prctile(irf_store_chol(:,:,:),97.5,3);
    irfs_lbb_chol(:,:)   = prctile(irf_store_chol(:,:,:),16,3);
    irfs_ubb_chol(:,:)   = prctile(irf_store_chol(:,:,:),84,3);

    irf_rr      =zeros(1,T);
    irf_rr_lb   =zeros(1,T);
    irf_rr_ub   =zeros(1,T);

    for jj=1:1:size(irf_rr,2)-12
        irf_rr(1,jj)    = irfs_chol(jj,index_i)     - (irfs_chol(jj+12,index_CPI)-irfs_chol(jj,index_CPI)) ;
        irf_rr_lb(1,jj) = irfs_lb_chol(jj,index_i)  - (irfs_lb_chol(jj+12,index_CPI)-irfs_lb_chol(jj,index_CPI)) ;
        irf_rr_ub(1,jj) = irfs_ub_chol(jj,index_i)  - (irfs_ub_chol(jj+12,index_CPI)-irfs_ub_chol(jj,index_CPI)) ;
    end

    Lambda=1/(1+0.04/12);
    qf      = zeros(size(irfs,2),1);

    for k=0:size(irfs,2)-110

        J=99;
        Funda=zeros(J+1,1);

        for j=0:J
            Funda(j+1,1)=(Lambda^j)*((1-Lambda)*irfs_chol(k+j+2,index_Div)-irf_rr(1,k+j+1)/12);
        end

        qf(k+1,1) = sum(Funda);

    end
    
    cd ..
    
    figure(1)
    for mm=1:size(irfs_chol,2)
        subplot(3,2,mm)
        plot(irfs_chol(:,mm),'Linewidth',3,'Color',rgb(irf_color))
        hold on
        plot(irfs_ub_chol(:,mm),'Linewidth',1,'Color',rgb(irf_color))
        plot(irfs_lb_chol(:,mm),'Linewidth',1,'Color',rgb(irf_color))
        plot(irfs_ubb_chol(:,mm),'Linewidth',1,'Color',rgb(irf_color))
        plot(irfs_lbb_chol(:,mm),'Linewidth',1,'Color',rgb(irf_color))
        ciplot(irfs_lb_chol(:,mm),irfs_ub_chol(:,mm), band_color);
        ciplot(irfs_lbb_chol(:,mm),irfs_ubb_chol(:,mm), band_color);
        if mm == 4
            plot(irf_rr(1,:),'Linewidth',3,'Linestyle','-.','Color',rgb(dash_color))
        elseif mm == 5
            plot(qf,'Linewidth',3,'Linestyle','-.','Color',rgb(dash_color))
        end
        plot(zeros(size(irfs(1,:),2),1),'k')
        axis tight
        xlim([1 IRF_T])
        ylim([y_lim(mm,1) y_lim(mm,2)])
        t=title(vars_chol_labels(1,mm));
        set(gca,'Box','on','FontSize', 16)
        ylabel('Percent','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
        if mm == 5 || mm == 6
            xlabel('Months','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
        end
        set(t,'fontsize',24,'FontName','Palatino','Interpreter','Latex')
        grid on
        set(gca,'GridLineStyle','--')
        ax = gca;
        ax.GridAlpha = .5;
    end
    h = gcf;
    set(h, 'PaperPositionMode','auto')
    set(h,'PaperOrientation','landscape')
    set(h,'Position',[0 0 1500 1000])
    tightfig
    print(h,'-depsc','..\Output\Figure_9')
end

%% Figure: Robustness with Excess Bond Premium

if choose_figure == 11
    
    figure(1)
    for mm=1:6
        subplot(3,2,mm)
        plot(irf_med(mm,:),'Linewidth',3,'Color',rgb(irf_color))
        hold on
        plot(irf_lb(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_ub(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_lbb(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        plot(irf_ubb(mm,:),'Linewidth',1,'Color',rgb(irf_color))
        ciplot(irf_lb(mm,:),irf_ub(mm,:), band_color);
        ciplot(irf_lbb(mm,:),irf_ubb(mm,:), band_color);
        if mm == 1
            plot(irfs(n_vars+1,:),'Linewidth',3,'Linestyle','-.','Color',rgb(dash_color))
        elseif mm == 2
            plot(qf,'Linewidth',3,'Linestyle','-.','Color',rgb(dash_color))
        end
        plot(zeros(size(irfs(mm,:),2),1),'k')
        axis tight
        xlim([1 IRF_T])
        ylim([y_lim(mm,1) y_lim(mm,2)])
        t=title(labels(1,mm));
        set(gca,'Box','on','FontSize', 16)
        ylabel('Percent','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
        if mm == 5 || mm == 6
            xlabel('Months','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
        end
        set(t,'fontsize',24,'FontName','Palatino','Interpreter','Latex')
        grid on
        set(gca,'GridLineStyle','--')
        ax = gca;
        ax.GridAlpha = .5;
    end
    h = gcf;
    set(h, 'PaperPositionMode','auto')
    set(h,'PaperOrientation','landscape')
    set(h,'Position',[0 0 1500 1000])
    tightfig
    print(h,'-depsc','..\Output\Figure_11_1')

    figure(2)
    subplot(3,2,1)
    plot(irf_med(7,:),'Linewidth',3,'Color',rgb(irf_color))
    hold on
    plot(irf_lb(7,:),'Linewidth',1,'Color',rgb(irf_color))
    plot(irf_ub(7,:),'Linewidth',1,'Color',rgb(irf_color))
    plot(irf_lbb(7,:),'Linewidth',1,'Color',rgb(irf_color))
    plot(irf_ubb(7,:),'Linewidth',1,'Color',rgb(irf_color))
    ciplot(irf_lb(7,:),irf_ub(7,:), band_color);
    ciplot(irf_lbb(7,:),irf_ubb(7,:), band_color);
    plot(zeros(size(irfs(1,:),2),1),'k')
    axis tight
    xlim([1 IRF_T])
    ylim([y_lim(7,1) y_lim(7,2)])
    t= title(labels(1,7));
    set(gca,'Box','on','FontSize', 16)
    ylabel('Percent','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
    xlabel('Months','Fontsize',20,'FontName','Palatino','Interpreter','Latex')
    set(t,'fontsize',24,'FontName','Palatino','Interpreter','Latex')
    grid on
    set(gca,'GridLineStyle','--')
    ax = gca;
    ax.GridAlpha = .5;
    h = gcf;
    set(h, 'PaperPositionMode','auto')
    set(h,'PaperOrientation','landscape')
    set(h,'Position',[0 0 1500 1000])
    tightfig
    print(h,'-depsc','..\Output\Figure_11_2')

end
